* D:\E\replications\BJPS2020\resultsets\figures_1_4_5.do
	// was called variability_december2017.do

* cd D:\E\replications\BJPS2020\
set more off
clear all
cap log close
log using ./resultsets/figures_1_4_5.txt, text replace
display "$S_TIME  $S_DATE"
about
ado dir


**********
* Figures 1, 4, and appendix-Table-1
**********

* Figure 1: Prevalence of large income losses
	* Appendix-Figure 1: Prevalence of large income losses (50%)
* Figure 4: Risk reduction (pre- and post-2008)
	* Appendix-Figure 3: Risk reduction (pre- and post-2008); 50% income losses

* Prep for figures 1 & 4 (and Appendix-Figure 1 & 3)  
clear
foreach j in 25 10 50 {
	use ./resultsets/resultset_drops.dta, clear
	keep if ages=="25-60"
	xtset ccode year

	gen year_di=year if ESI`j'_di_e<.
	gen year_mi=year if ESI`j'_mi_e<.
	gen year_red=year if red2_`j'<.

	*local j 25 // 10 25 50
	* MI drops
	preserve
		encode dataset, gen(d)
		statsby _b _se e(df_r), by(ccode) clear : reg ESI`j'_mi_e year i.d
		isid ccode
		gen t=_b_year/_se_year
		gen df=_eq2_stat_1
		gen p=2*ttail(df, abs(t))
		gen trend="up" if t>0 & t<. & p<=0.05
		replace trend="down" if t<0 & p<=0.05
		replace trend="none" if p>0.05 & p<.
		list ccode trend
		keep ccode trend
		rename trend trend_mi
		tempfile t
		save `t'
	restore

	* Di drops
	preserve
		encode dataset, gen(d)
		statsby _b _se e(df_r), by(ccode) clear : reg ESI`j'_di_e year i.d
		isid ccode
		gen t=_b_year/_se_year
		gen df=_eq2_stat_1
		gen p=2*ttail(df, abs(t))
		gen trend="up" if t>0 & t<. & p<=0.05
		replace trend="down" if t<0 & p<=0.05
		replace trend="none" if p>0.05 & p<.
		list ccode trend p
		keep ccode trend
		rename trend trend_di
		tempfile t2
		save `t2'
	restore

	* Risk reduction
	preserve
		encode dataset, gen(d)
		statsby _b _se e(df_r), by(ccode) clear : reg red2_`j' year i.d
		isid ccode
		gen t=_b_year/_se_year
		gen df=_eq2_stat_1
		gen p=2*ttail(df, abs(t))
		gen trend="up" if t>0 & t<. & p<=0.05
		replace trend="down" if t<0 & p<=0.05
		replace trend="none" if p>0.05 & p<.
		list ccode trend p t
		keep ccode trend
		rename trend trend_riskreduction
		tempfile t3
		save `t3'
	restore

	set more off	
	foreach k of varlist year year_di year_mi year_red red2_`j' ESI??_mi_e ESI??_di_e { // d`j'_mi_e d`j'_di_e { // gini_mi_e gini_di_e {
		foreach v in min max my p50 {
			clonevar `v'_`k'=`k'
		}
		foreach v in pre {
			clonevar `v'_`k'=`k' if year<=2008
		}
		foreach v in post {
			clonevar `v'_`k'=`k' if year>2008 & year<.
		}
	}

	* collapse (min) min_* (max) max_* (mean) pre_* post_* my_*  rank_* (p50) p50_*, by(ccode iso3c)
	collapse (min) min_* (max) max_* (mean) pre_* post_* my_*  (p50) p50_*, by(ccode iso3c)
	merge 1:1 ccode using `t'
	drop _merge
	merge 1:1 ccode using `t2'
	cap drop _merge
	merge 1:1 ccode using `t3'
	drop _merge

	cap drop group
	egen group_MI=group(my_ESI`j'_mi_e) 
	egen group_DI=group(my_ESI`j'_di_e) 
	cap drop group_riskreduction
	egen group_riskreduction=group(my_red2_`j')
	
	cap drop cyr
	gen str cyr=iso3c+" ("+substr(string(min_year_mi),-2,2)+"-"+substr(string(max_year_mi),-2,2)+")"
	labmask group_MI, val(cyr)
	labmask group_DI, val(cyr)
	labmask group_riskreduction, val(cyr)

	cap drop group1
	gen group1=group_MI+.125
	cap drop group2
	gen group2=group_MI-.125

	* Figure 1
	* local j 25
	sum ccode
	local m `r(N)'
	* local j 25
	sum max_ESI`j'_mi_e
	local xl=round(`r(max)',.1)
	di "`xl'"
	twoway  (scatter group_MI my_ESI`j'_mi_e, sort ms(none)) ///
			(scatter group1   my_ESI`j'_mi_e, sort ms(o)) ///
			(scatter group2   my_ESI`j'_di_e, sort ms(oh)) ///
			(pcarrow group1 min_ESI`j'_mi_e group1  max_ESI`j'_mi_e if trend_mi=="up", msize(medlarge)) ///
			(pcarrow group1 max_ESI`j'_mi_e group1  min_ESI`j'_mi_e if trend_mi=="down", msize(medlarge)) ///		
			(pcarrow group1 min_ESI`j'_mi_e group1  max_ESI`j'_mi_e if trend_mi=="none", mstyle(none) ) ///
			(pcarrow group2 min_ESI`j'_di_e group2  max_ESI`j'_di_e if trend_di=="up"  , lp(dash) msize(medlarge)) ///
			(pcarrow group2 max_ESI`j'_di_e group2  min_ESI`j'_di_e if trend_di=="down", lp(dash) msize(medlarge)) ///		
			(pcarrow group2 min_ESI`j'_di_e group2  max_ESI`j'_di_e if trend_di=="none", lp(dash) mstyle(none)) ///
			, scheme(rbn1mono) /// ylabel(1(1)19) ///
			yla(1/`m', ang(h) valuelabel) /// notick labsize(*0.75)) ///
			ytitle("") xtitle(" " "Incidence of `j'%+ arc income losses") ///
			/// xlabel(, angle(0)) ///
			xlabel(0(.05)`xl', angle(0) format(%3.2f)) ///
			/// legend(holes(4) label(1 "") label(2 "MI") label(3 "min-max MI") label(4 "DI") label(5 "min-max DI") r(2) size(small)) ///
			legend(r(2) order(2 "Market income (mean)"  6 "min-max (w/ arrow for trend)" 3 "Disposable income (mean)" 9 "min-max (w/ arrow for trend)")) ///
			legend(region(lwidth(none))) ///
			ysize(4.5) xsize(3) 

	graph save   ./resultsets/figures/figure1_`j'.gph, replace
	*graph export ./resultsets/figures/figure1_`j'.emf, replace
	graph export ./resultsets/figures/figure1_`j'.pdf, replace
	
	* Figure 4
	foreach v of varlist my_red2_?? min_red2_?? max_red2_?? pre_red2_?? post_red2_?? {
		replace `v'=`v'*100
	}

	cap drop cyr
	gen str cyr=iso3c+" ("+substr(string(min_year_mi),-2,2)+"-"+substr(string(max_year_mi),-2,2)+")"
	labmask group_MI, val(cyr)
	labmask group_DI, val(cyr)
	labmask group_riskreduction, val(cyr)
		
	* local j 25
	* local j 50
	sum ccode
	local m `r(N)'
	* local j 50
	sum max_red2_`j'
	local xl=round(`r(max)',.1)
	di "`xl'"
	twoway  (scatter group_riskreduction   my_red2_`j', sort ms(none)) ///
			(scatter group_riskreduction   pre_red2_`j', sort ms(X)) ///
			(pcarrow group_riskreduction min_red2_`j' group_riskreduction  max_red2_`j' if trend_riskreduction=="up"  , lp(solid) msize(medlarge)) ///
			(pcarrow group_riskreduction max_red2_`j' group_riskreduction  min_red2_`j' if trend_riskreduction=="down", lp(solid) msize(medlarge)) ///		
			(pcarrow group_riskreduction min_red2_`j' group_riskreduction  max_red2_`j' if trend_riskreduction=="none", lp(solid) mstyle(none)) ///
			(scatter group_riskreduction   post_red2_`j', sort ms(Oh)) ///
			, scheme(rbn1mono) /// ylabel(1(1)19) ///
			yla(1/`m', ang(h) valuelabel) /// notick labsize(*0.75)) ///
			ytitle("") xtitle(" " "Risk-reduction (`j'+% arc income losses)") ///
			/// xlabel(, angle(0) format(%2.0f)) ///
			xlabel(0(5)`xl', angle(0) format(%2.0f)) ///
			/// legend(holes(4) label(1 "") label(2 "MI") label(3 "min-max MI") label(4 "DI") label(5 "min-max DI") r(2) size(small)) ///
			legend(r(4) order(2 "Risk-reduction (mean, pre-2008)" 6 "Risk-reduction (mean, post-2008)" 5 "Range (min-max), no trend" 3 "Range (min-max), w/ trend arrow")) ///
			legend(region(lwidth(none))) ///
			ysize(4.5) xsize(3) 

	graph save   ./resultsets/figures/figure4_`j'.gph, replace
	*graph export ./resultsets/figures/figure4_`j'.emf, replace
	graph export ./resultsets/figures/figure4_`j'.pdf, replace
}

****************
* Figure 5
****************

* Figure 5: Social spending, as well as benefit generosity, and risk reduction
	* Appendix-Figure 4: Social spending, as well as benefit generosity, and risk reduction (50%)

use ./resultsets/resultset_drops.dta, clear
keep if ages=="25-60"
xtset ccode year

	bys ccode (year): gen first=(_n==1)
	label var ls_totgen "Combined Generosity Index"
	label var ls_uegen  "Unemployment Generosity Index"
	label var ls_skgen  "Sickness Generosity Index"
	label var ls_pgen   "Pension Generosity Index"
	label var M_ls_totgen "Combined Generosity Index"
	label var M_ls_uegen  "Unemployment Generosity Index"
	label var M_ls_skgen  "Sickness Generosity Index"
	label var M_ls_pgen   "Pension Generosity Index"
	label var socx_all    "Total public social expenditure (% GDP)"
	label var M_socx_all  "Total public social expenditure (% GDP)"

	format *ls_* %2.0f
	format *red2* %2.1f
	
	foreach j in 10 25 50 {
		* local j 25
		label var red2_`j' "Risk-reduction (`j'+% arc income losses)"
		label var M_red2_`j' "Risk-reduction (`j'+% arc income losses)"
		local dv2 red2_`j'
		local dv M_`dv2'
		foreach rhs2 of varlist socx_all ls_totgen ls_uegen ls_skgen  { 
			local rhs M_`rhs2'
			reg `dv' `rhs' if first
			local C=round(_b[`rhs'],0.001)
			local c=substr(string(`C'),1,5)
			local SE=round(_se[`rhs'],0.001)
			local se=substr(string(`SE'),1,5)
			local n=substr(string(`e(N)'),1,4)
			local T=round(`=_b[`rhs']/_se[`rhs']',0.1)
			local t=substr(string(`T'),1,4)
			local R=round(`e(r2_a)',0.01)
			local r=substr(string(`R'),1,4)
			di "`c', `se', `t', `r', `n'"
			local d
			levelsof ccode, local(country)
				foreach cnty of local country {
				local d "`d' (lfit `dv2' `rhs2' if ccode==`cnty', lp(solid) lw(thin) lc(gs9))"
			}
			tw  ///
				(scatter `dv2' `rhs2', ms(oh) mcolor(gs13) msize(vsmall)) /// mlabcolor(gs9) mlab(iso3c) mlabsize(tiny) mlabpos(0)
				(lfit `dv2' `rhs2', lp(dash) lc(gs13)) ///
				(lfit `dv' `rhs' if first) ///
				(scatter `dv' `rhs' if first, mlab(iso3c) mlabpos(0) ms(none)) ///
				, xtitle("`: var label `rhs''") ///
				ytitle("`: var label `dv''") ///
				scheme(rbn1mono) ///
				xlabel(, angle(0)) ylabel(, angle(0)) ///
				legend(off) ///
				note("B/w effects regression (solid line):" "Coef: `c', SE: `se', T: `t', R2: `r', N: `n'", size(vsmall))

			graph save ./resultsets/figures/tmp/tmpFig5_`dv'_`rhs'_`j'.gph, replace	
			more
		}

	
	* Figure 5
	graph combine ///
		./resultsets/figures/tmp/tmpFig5_M_red2_`j'_M_socx_all_`j'.gph ///
		./resultsets/figures/tmp/tmpFig5_M_red2_`j'_M_ls_totgen_`j'.gph ///
		./resultsets/figures/tmp/tmpFig5_M_red2_`j'_M_ls_uegen_`j'.gph ///
		./resultsets/figures/tmp/tmpFig5_M_red2_`j'_M_ls_skgen_`j'.gph ///
		, ysize(6) xsize(5) scheme(rbn1mono)
	
	graph save   ./resultsets/figures/figure5_`j'.gph, replace
	*graph export ./resultsets/figures/figure5_`j'.emf, replace
	graph export ./resultsets/figures/figure5_`j'.pdf, replace
	}
	
******************

* Appendix-Table 1: Correlation between of key variables with different cut-off points (10/25/50%)
use ./resultsets/resultset_drops.dta, clear
keep if ages=="25-60"
xtset ccode year

	pwcorr ESI10_mi_e ESI25_mi_e ESI50_mi_e ESI10_di_e ESI25_di_e ESI50_di_e RR10_2 RR25_2 RR50_2, star(0.05)
	* This is Appendix-Table 1:
	pwcorr ESI10_mi_e ESI25_mi_e ESI50_mi_e ESI10_di_e ESI25_di_e ESI50_di_e RR10_2 RR25_2 RR50_2


display "$S_TIME  $S_DATE"
cap log close

